Preliminaries
Loading exec control data
Create folds
Entropy PCA and scaling of variables
Join LHQ and AXCPT
Jobs data
Experiment version
Data plots
Set loo fold
Traditional AXCPT analysis

Preliminaries

Loading exec control data

exec_lhq <- as.data.table(fread("subjectdata/axcpt_lhq_data_forpub.csv"))
axcpt <- as.data.table(fread("rawdata/axcpt_trial_data.csv", stringsAsFactors = F, 
    sep = ","))

Create folds

k = 10
k_loo = length(unique(exec_lhq$participant))

exec_lhq <- fold(exec_lhq, k = k, id_col = "participant") %>% ungroup()
exec_lhq <- fold(exec_lhq, k = k_loo, id_col = "participant") %>% ungroup()

## Warning in fold(exec_lhq, k = k_loo, id_col = "participant"): Found 1 existing
## fold column. It will NOT be replaced. Change 'handle_existing_fold_cols' to
## 'remove' if you want to replace it.

exec_lhq <- exec_lhq %>% rename(.folds.10 = .folds_1, .folds.loo = .folds_2)

Entropy PCA and scaling of variables

exec_lhq_ppd <- preprocess_steps(exec_lhq)

Join LHQ and AXCPT

# Merge behavioral and subject data (restricting to only those in subject data
# frame)
axcpt_ppd <- left_join(axcpt, exec_lhq_ppd, by = "participant", suffix = c(".x", 
    ""))
axcpt_ppd <- axcpt_preprocess(axcpt_ppd)

sd.rt <- axcpt_ppd %>% group_by(participant) %>% summarise(meanRT = mean(RT_correct, 
    na.rm = T)) %>% ungroup()

sd.acc <- axcpt_ppd %>% group_by(participant) %>% summarise(meanAcc = mean(accuracy_target))

sd.rt <- sd.rt %>% mutate(gmean = mean(meanRT), gsd = sd(meanRT), lower = gmean - 
    2.5 * gsd, upper = gmean + 2.5 * gsd, outlier = ifelse(meanRT >= lower & meanRT <= 
    upper, "fine", "outlier"))

exclude <- sd.rt$participant[sd.rt$outlier == "outlier"]

axcpt_ppd <- axcpt_ppd[!(axcpt_ppd$participant %in% exclude), ]

exec_lhq <- exec_lhq[!(exec_lhq$participant %in% exclude), ]
exec_lhq_ppd <- exec_lhq_ppd[!(exec_lhq_ppd$participant %in% exclude), ]

exec_lhq_ppd <- fold(exec_lhq_ppd, k = k_loo - length(exclude), id_col = "participant", 
    handle_existing_fold_cols = "remove") %>% ungroup()

k_loo = k_loo - length(exclude)

exec_lhq_ppd <- exec_lhq_ppd %>% rename(.folds.loo = .folds)

newfolds <- exec_lhq_ppd %>% select(participant, .folds.loo)

axcpt_ppd <- axcpt_ppd %>% select(-.folds.loo)
axcpt_ppd <- axcpt_ppd %>% left_join(newfolds)

## Joining, by = "participant"

Jobs data

axcpt_ppd$occupation <- factor(axcpt_ppd$occupation)
contrasts(axcpt_ppd$occupation) <- contr.sum(4)

Experiment version

axcpt_ppd$ExperimentName <- factor(axcpt_ppd$ExperimentName)
contrasts(axcpt_ppd$ExperimentName) <- contr.sum(3)

Data plots

Summary stats table

data.agg <- exec_lhq %>% select(participant, age, aoa, yearsBilingual, L1_current_exposure, 
    L2_current_exposure, L3_current_exposure, ReadEntropy, SpeakEntropy, HomeEntropy, 
    WorkEntropy, SocialEntropy, meanAccent, meanL2ability)

stargazer(data.agg, type = "html", out = "tables/summary_stats.html")

## 
## <table style="text-align:center"><tr><td colspan="8" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Pctl(25)</td><td>Pctl(75)</td><td>Max</td></tr>
## <tr><td colspan="8" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">participant</td><td>455</td><td>758.965</td><td>323.781</td><td>212</td><td>453.5</td><td>997</td><td>1,264</td></tr>
## <tr><td style="text-align:left">age</td><td>455</td><td>22.754</td><td>3.633</td><td>18</td><td>20</td><td>24</td><td>35</td></tr>
## <tr><td style="text-align:left">aoa</td><td>455</td><td>6.752</td><td>4.334</td><td>0</td><td>4</td><td>10</td><td>28</td></tr>
## <tr><td style="text-align:left">yearsBilingual</td><td>455</td><td>16.002</td><td>5.340</td><td>1</td><td>13</td><td>19.6</td><td>34</td></tr>
## <tr><td style="text-align:left">L1_current_exposure</td><td>455</td><td>63.690</td><td>20.576</td><td>5</td><td>50</td><td>80</td><td>98</td></tr>
## <tr><td style="text-align:left">L2_current_exposure</td><td>455</td><td>33.890</td><td>20.070</td><td>2</td><td>20</td><td>50</td><td>95</td></tr>
## <tr><td style="text-align:left">L3_current_exposure</td><td>455</td><td>2.420</td><td>5.648</td><td>0</td><td>0</td><td>1</td><td>40</td></tr>
## <tr><td style="text-align:left">ReadEntropy</td><td>455</td><td>0.595</td><td>0.410</td><td>0.000</td><td>0.286</td><td>0.971</td><td>1.571</td></tr>
## <tr><td style="text-align:left">SpeakEntropy</td><td>455</td><td>0.695</td><td>0.417</td><td>0.000</td><td>0.469</td><td>1.000</td><td>1.585</td></tr>
## <tr><td style="text-align:left">HomeEntropy</td><td>455</td><td>0.602</td><td>0.456</td><td>0.000</td><td>0.000</td><td>0.963</td><td>1.585</td></tr>
## <tr><td style="text-align:left">WorkEntropy</td><td>455</td><td>0.756</td><td>0.367</td><td>0</td><td>0.6</td><td>1.0</td><td>2</td></tr>
## <tr><td style="text-align:left">SocialEntropy</td><td>455</td><td>0.944</td><td>0.276</td><td>0.000</td><td>0.811</td><td>1.000</td><td>1.585</td></tr>
## <tr><td style="text-align:left">meanAccent</td><td>455</td><td>3.625</td><td>1.787</td><td>1.000</td><td>2.000</td><td>5.000</td><td>7.000</td></tr>
## <tr><td style="text-align:left">meanL2ability</td><td>455</td><td>7.319</td><td>1.847</td><td>1.000</td><td>6.300</td><td>8.800</td><td>10.000</td></tr>
## <tr><td colspan="8" style="border-bottom: 1px solid black"></td></tr></table>

Group level

sd.rt <- sem(axcpt_ppd, RT_correct, participant, Condition)
sd.acc <- sem(axcpt_ppd, accuracy_target, participant, Condition)

sd.rt$Condition <- factor(sd.rt$Condition, levels = c("AY", "BX", "BY", "AX"))
sd.acc$Condition <- factor(sd.acc$Condition, levels = c("AY", "BX", "BY", "AX"))

p1 <- ggplot(sd.rt, aes(x = Condition, y = mean_RT_correct, ymin = lower, ymax = upper, 
    fill = Condition)) + geom_col() + geom_errorbar(width = 0.2) + ylab("Mean RT (in ms)") + 
    xlab("Trial type") + labs(fill = "Trial type") + theme_minimal(base_size = 12) + 
    scale_fill_brewer(palette = "Set1")


p2 <- ggplot(sd.acc, aes(x = Condition, y = mean_accuracy_target, ymin = lower, ymax = upper, 
    fill = Condition)) + geom_col() + geom_errorbar(width = 0.2) + ylab("Mean accuracy") + 
    xlab("Trial type") + labs(fill = "Trial type") + theme_minimal(base_size = 12) + 
    scale_fill_brewer(palette = "Set1")

(plot <- ggarrange(p1, p2, ncol = 2, labels = c("A", "B"), common.legend = T, legend = "right", 
    align = "v"))
ggsave(plot = plot, file = "figures/mean_acc_rt.png", width = 12, height = 6, units = "cm")

RT: individual differences

rt.sd <- axcpt_ppd %>% group_by(participant, L2_current_exposure, aoa, PCA_General, 
    PCA_Work, Condition) %>% summarise(meanRT = mean(RT_correct, na.rm = T), sdRT = sd(RT_correct, 
    na.rm = T), N = n(), semRT = sdRT/sqrt(N))


rt.sd$Condition <- factor(rt.sd$Condition, levels = c("AY", "BX", "BY", "AX"))

p1 <- ggplot(rt.sd, aes(x = aoa, y = meanRT, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Reaction time (in ms)") + xlab("L2 AoA") + theme_minimal(base_size = 12) + 
    scale_colour_brewer(palette = "Set1")

p2 <- ggplot(rt.sd, aes(x = L2_current_exposure, y = meanRT, colour = Condition)) + 
    geom_point(alpha = 0.2) + geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + 
    xlab("L2 current exposure") + theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p3 <- ggplot(rt.sd, aes(x = PCA_General, y = meanRT, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Reaction time (in ms)") + xlab("PC: General entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p4 <- ggplot(rt.sd, aes(x = PCA_Work, y = meanRT, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Reaction time (in ms)") + xlab("PC: Work entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")


(plot <- ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2, labels = c("A", "B", "C", 
    "D"), common.legend = T, legend = "right", align = "hv"))
ggsave(plot = plot, file = "figures/rt_by_id.png", width = 17.8, height = 17.5, units = "cm")

RT: individual differences: raw

raw <- axcpt_ppd[]
raw$Condition <- factor(raw$Condition, levels = c("AY", "BX", "BY", "AX"))



p1 <- ggplot(raw, aes(x = aoa, y = RT_correct, colour = Condition)) + geom_point(alpha = 0.05) + 
    geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + xlab("L2 AoA") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p2 <- ggplot(raw, aes(x = L2_current_exposure, y = RT_correct, colour = Condition)) + 
    geom_point(alpha = 0.05) + geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + 
    xlab("L2 current exposure") + theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p3 <- ggplot(raw, aes(x = PCA_General, y = RT_correct, colour = Condition)) + geom_point(alpha = 0.05) + 
    geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + xlab("PC: General entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p4 <- ggplot(raw, aes(x = PCA_Work, y = RT_correct, colour = Condition)) + geom_point(alpha = 0.05) + 
    geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + xlab("PC: Work entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")


(plot <- ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2, labels = c("A", "B", "C", 
    "D"), common.legend = T, legend = "right", align = "hv"))

## Warning: Removed 5051 rows containing non-finite values (stat_smooth).

## Warning: Removed 5051 rows containing missing values (geom_point).

## Warning: Removed 5051 rows containing non-finite values (stat_smooth).

## Warning: Removed 5051 rows containing missing values (geom_point).

## Warning: Removed 5051 rows containing non-finite values (stat_smooth).

## Warning: Removed 5051 rows containing missing values (geom_point).

## Warning: Removed 5051 rows containing non-finite values (stat_smooth).

## Warning: Removed 5051 rows containing missing values (geom_point).

## Warning: Removed 5051 rows containing non-finite values (stat_smooth).

## Warning: Removed 5051 rows containing missing values (geom_point).

ggsave(plot = plot, file = "figures/rt_by_id_raw.png", width = 17.8, height = 17.5, 
    units = "cm")

ACC: individual differences

acc.sd <- axcpt_ppd %>% group_by(participant, L2_current_exposure, aoa, PCA_General, 
    PCA_Work, Condition) %>% summarise(meanACC = mean(accuracy_target, na.rm = T), 
    sdACC = sd(accuracy_target, na.rm = T), N = n(), semACC = sdACC/sqrt(N))

p1 <- ggplot(acc.sd, aes(x = aoa, y = meanACC, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Mean accuracy") + xlab("L2 AoA") + theme_minimal(base_size = 12) + 
    scale_colour_brewer(palette = "Set1")

p2 <- ggplot(acc.sd, aes(x = L2_current_exposure, y = meanACC, colour = Condition)) + 
    geom_point(alpha = 0.2) + geom_smooth(method = "lm") + ylab("Mean accuracy") + 
    xlab("L2 current exposure") + theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p3 <- ggplot(acc.sd, aes(x = PCA_General, y = meanACC, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Mean accuracy") + xlab("PC: General entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p4 <- ggplot(acc.sd, aes(x = PCA_Work, y = meanACC, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Mean accuracy") + xlab("PC: Work entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")


(plot <- ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2, labels = c("A", "B", "C", 
    "D"), common.legend = T, legend = "right", align = "hv"))
ggsave(plot = plot, file = "figures/acc_by_id.png", width = 17.8, height = 17.5, 
    units = "cm")

Set loo fold

axcpt_ppd$.folds <- axcpt_ppd$.folds.loo
axcpt_ppd <- axcpt_ppd[axcpt_ppd$Condition != "AX", ]

Traditional AXCPT analysis

mod.rt.0 <- "RT_correct  ~  Condition + (1|participant)"

mod.rt.1 <- "RT_correct  ~  Condition*aoa + Condition*L2_exposure +  (1|participant)"

mod.rt.1.1 <- "RT_correct  ~  Condition*PCA_Work + Condition*PCA_General +  (1|participant)"

mod.rt.2 <- "RT_correct  ~  Condition*aoa + Condition*L2_exposure + Condition*PCA_Work + Condition*PCA_General +  (1|participant)"

mod.rt.3 <- "RT_correct  ~  aoa * (Condition*L2_exposure + Condition*PCA_Work + Condition*PCA_General) +  (1|participant)"



mod.acc.0 <- "acc        ~   Condition + (1|participant)"

mod.acc.0.1 <- "acc        ~   Condition + L2_exposure*aoa + (1|participant)"

mod.acc.1 <- "acc        ~   Condition*aoa + Condition*L2_exposure + (1|participant)"

mod.acc.1.1 <- "acc        ~   Condition*PCA_Work + Condition*PCA_General + (1|participant)"

mod.acc.2 <- "acc        ~   Condition*aoa + Condition*L2_exposure + Condition*PCA_Work + Condition*PCA_General + (1|participant)"

mod.acc.3 <- "acc        ~   aoa * (Condition*L2_exposure + Condition*PCA_Work + Condition*PCA_General) + (1|participant)"

RT models

# axcpt_ppd <- axcpt_ppd[!(axcpt_ppd$participant %in% exclude),]
axcpt_ppd$Condition <- relevel(axcpt_ppd$Condition, ref = "BX")
all.rt.0 <- lmer(mod.rt.0, axcpt_ppd, REML = F)
all.rt.1 <- lmer(mod.rt.1, axcpt_ppd, REML = F)
all.rt.1.1 <- lmer(mod.rt.1.1, axcpt_ppd, REML = F)
all.rt.2 <- lmer(mod.rt.2, axcpt_ppd, REML = F)
all.rt.3 <- lmer(mod.rt.3, axcpt_ppd, REML = F)

anova(all.rt.0, all.rt.1, all.rt.2, all.rt.3)

## Data: axcpt_ppd
## Models:
## all.rt.0: mod.rt.0
## all.rt.1: mod.rt.1
## all.rt.2: mod.rt.2
## all.rt.3: mod.rt.3
##          Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)  
## all.rt.0  5 254753 254793 -127372   254743                            
## all.rt.1 11 254754 254841 -127366   254732 11.0123      6    0.08800 .
## all.rt.2 17 254752 254886 -127359   254718 13.8777      6    0.03103 *
## all.rt.3 26 254764 254969 -127356   254712  6.0874      9    0.73114  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kable(anova(all.rt.0, all.rt.1, all.rt.2, all.rt.3))

	Df	AIC	BIC	logLik	deviance	Chisq	Chi Df	Pr(>Chisq)
all.rt.0	5	254753.1	254792.6	-127371.6	254743.1	NA	NA	NA
all.rt.1	11	254754.1	254840.9	-127366.1	254732.1	11.012288	6	0.0879974
all.rt.2	17	254752.2	254886.4	-127359.1	254718.2	13.877692	6	0.0310327
all.rt.3	26	254764.2	254969.4	-127356.1	254712.2	6.087412	9	0.7311388

kable(anova(all.rt.0, all.rt.1.1, all.rt.2, all.rt.3))

	Df	AIC	BIC	logLik	deviance	Chisq	Chi Df	Pr(>Chisq)
all.rt.0	5	254753.1	254792.6	-127371.6	254743.1	NA	NA	NA
all.rt.1.1	11	254746.0	254832.9	-127362.0	254724.0	19.095428	6	0.0040055
all.rt.2	17	254752.2	254886.4	-127359.1	254718.2	5.794552	6	0.4465938
all.rt.3	26	254764.2	254969.4	-127356.1	254712.2	6.087412	9	0.7311388

AIC(all.rt.0, all.rt.1, all.rt.1.1, all.rt.2, all.rt.3)

##            df      AIC
## all.rt.0    5 254753.1
## all.rt.1   11 254754.1
## all.rt.1.1 11 254746.0
## all.rt.2   17 254752.2
## all.rt.3   26 254764.2

BIC(all.rt.0, all.rt.1, all.rt.1.1, all.rt.2, all.rt.3)

##            df      BIC
## all.rt.0    5 254792.6
## all.rt.1   11 254840.9
## all.rt.1.1 11 254832.9
## all.rt.2   17 254886.4
## all.rt.3   26 254969.4

So entropy explains more then L2 exp and aoa but not vice versa

summ(all.rt.1.1, confint = T, digits = 3)

Observations	19796
Dependent variable	RT_correct
Type	Mixed effects linear regression

AIC	254746.033
BIC	254832.859
Pseudo-R² (fixed effects)	0.136
Pseudo-R² (total)	0.315

Fixed Effects
	Est.	2.5%	97.5%	t val.	d.f.	p
(Intercept)	384.985	377.229	392.742	97.283	607.866	0.000
ConditionAY	113.831	108.715	118.947	43.608	19364.317	0.000
ConditionBY	-42.140	-47.041	-37.240	-16.854	19351.923	0.000
PCA_Work	9.540	1.337	17.743	2.279	608.076	0.023
PCA_General	-4.424	-12.618	3.770	-1.058	605.252	0.290
ConditionAY:PCA_Work	-5.449	-10.860	-0.037	-1.974	19361.250	0.048
ConditionBY:PCA_Work	-1.739	-6.924	3.446	-0.657	19351.850	0.511
ConditionAY:PCA_General	9.271	3.889	14.654	3.376	19363.649	0.001
ConditionBY:PCA_General	7.448	2.288	12.608	2.829	19352.413	0.005
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
participant	(Intercept)	74.973
Residual		146.380

Grouping Variables
Group	# groups	ICC
participant	455	0.208

summary(all.rt.1.1)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: mod.rt.1.1
##    Data: axcpt_ppd
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  254746.0  254832.9 -127362.0  254724.0     19785 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1265 -0.6608 -0.1892  0.4511  4.8437 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  participant (Intercept)  5621     74.97  
##  Residual                21427    146.38  
## Number of obs: 19796, groups:  participant, 455
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)               384.985      3.957   607.866  97.283  < 2e-16 ***
## ConditionAY               113.831      2.610 19364.317  43.608  < 2e-16 ***
## ConditionBY               -42.140      2.500 19351.923 -16.854  < 2e-16 ***
## PCA_Work                    9.540      4.185   608.076   2.279 0.022984 *  
## PCA_General                -4.424      4.181   605.252  -1.058 0.290400    
## ConditionAY:PCA_Work       -5.449      2.761 19361.250  -1.974 0.048452 *  
## ConditionBY:PCA_Work       -1.739      2.646 19351.850  -0.657 0.511042    
## ConditionAY:PCA_General     9.271      2.746 19363.649   3.376 0.000737 ***
## ConditionBY:PCA_General     7.448      2.633 19352.413   2.829 0.004675 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) CndtAY CndtBY PCA_Wr PCA_Gn CAY:PCA_W CBY:PCA_W CAY:PCA_G
## ConditionAY -0.313                                                          
## ConditionBY -0.326  0.493                                                   
## PCA_Work    -0.010  0.013  0.014                                            
## PCA_General  0.005 -0.007 -0.008 -0.331                                     
## CndAY:PCA_W  0.013 -0.041 -0.021 -0.313  0.103                              
## CndBY:PCA_W  0.014 -0.021 -0.042 -0.326  0.107  0.493                       
## CndAY:PCA_G -0.007  0.018  0.011  0.103 -0.311 -0.333    -0.162             
## CndBY:PCA_G -0.008  0.011  0.025  0.107 -0.324 -0.162    -0.331     0.493

anova(all.rt.1.1)

## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq  Mean Sq NumDF   DenDF   F value    Pr(>F)    
## Condition             82348357 41174178     2 19356.6 1921.5931 < 2.2e-16 ***
## PCA_Work                 72307    72307     1   453.9    3.3746  0.066864 .  
## PCA_General               1873     1873     1   453.0    0.0874  0.767652    
## Condition:PCA_Work       86262    43131     2 19355.2    2.0129  0.133627    
## Condition:PCA_General   282581   141291     2 19356.4    6.5940  0.001372 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ef <- as.data.frame(Effect(c("PCA_Work", "Condition"), all.rt.1.1, xlevels = list(PCA_Work = c(-2, 
    -1, 0, 1, 2)), confidence.level = 0.686))

f1a <- ggplot(ef, aes(x = PCA_Work, y = fit, ymin = lower, ymax = upper, group = Condition)) + 
    geom_line(aes(colour = Condition), size = 1) + geom_ribbon(fill = "grey", alpha = 0.3) + 
    coord_cartesian(ylim = c(300, 550)) + ylab("Model-estimated\nreaction time (in ms)") + 
    xlab("PC: Work Entropy") + scale_colour_brewer(palette = "Set1") + theme_bw(base_size = 12) + 
    labs(color = "Trial type") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

# ggsave('figures/axcpt_rt_work_condition.png', width=8.9, height=6, units =
# 'cm')

ef <- as.data.frame(Effect(c("PCA_General", "Condition"), all.rt.1.1, xlevels = list(PCA_General = c(-2, 
    -1, 0, 1, 2)), confidence.level = 0.686))

f1b <- ggplot(ef, aes(x = PCA_General, y = fit, ymin = lower, ymax = upper, group = Condition)) + 
    geom_line(aes(colour = Condition), size = 1) + geom_ribbon(fill = "grey", alpha = 0.3) + 
    coord_cartesian(ylim = c(300, 550)) + ylab("Model-estimates\nreaction time (in ms)") + 
    xlab("PC: General Entropy") + labs(color = "Trial type") + scale_colour_brewer(palette = "Set1") + 
    theme_bw(base_size = 12) + theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

# ggsave('figures/axcpt_rt_general_condition.png', width=8.9, height=6, units =
# 'cm')

(f1 <- ggarrange(f1b, f1a, labels = c("A", "B"), legend = "bottom", common.legend = T))

ggsave(plot = f1, "figures/traditional_rt_lme4_effects.png", width = 17.8, height = 8, 
    units = "cm")

RT models cross validated

cv.all.rt.0 <- crossvalidate(mod = mod.rt.0, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)

## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"

cv.all.rt.1 <- crossvalidate(mod = mod.rt.1, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)

## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"

cv.all.rt.1.1 <- crossvalidate(mod = mod.rt.1.1, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)

## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"

cv.all.rt.2 <- crossvalidate(mod = mod.rt.2, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)

## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"

cv.all.rt.3 <- crossvalidate(mod = mod.rt.3, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)

## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"

cv.rt.0.iter = cv.all.rt.0$performances
cv.rt.0 = cv.all.rt.0$performances_mean

cv.rt.1.iter = cv.all.rt.1$performances
cv.rt.1 = cv.all.rt.1$performances_mean

cv.rt.1.1.iter = cv.all.rt.1.1$performances
cv.rt.1.1 = cv.all.rt.1.1$performances_mean

cv.rt.2.iter = cv.all.rt.2$performances
cv.rt.2 = cv.all.rt.2$performances_mean

cv.rt.3.iter = cv.all.rt.3$performances
cv.rt.3 = cv.all.rt.3$performances_mean

xval_axcpt_rt_k_loo.iter <- rbind.data.frame(cv.rt.0.iter, cv.rt.1.iter, cv.rt.1.1.iter, 
    cv.rt.2.iter, cv.rt.3.iter)
write.csv(xval_axcpt_rt_k_loo.iter, "tables/xval_axcpt_rt_k_loo.iter.csv", row.names = F)

xval_axcpt_rt_k_loo.iter <- xval_axcpt_rt_k_loo.iter %>% group_by(mod) %>% mutate(meanRMSE_test = mean(RMSE_test)) %>% 
    ungroup() %>% mutate(mod = fct_reorder(mod, desc(meanRMSE_test)))

mean.xval_axcpt_rt_k_loo.iter <- xval_axcpt_rt_k_loo.iter %>% group_by(mod) %>% summarise(meanRMSE_train = mean(RMSE_train), 
    meanRMSE_test = mean(RMSE_test), meanr2_train = mean(r2_train), meanr2_test = mean(r2_test), 
    sd = sd(RMSE_test), N = n(), sem = sd/sqrt(N)) %>% ungroup() %>% mutate(mod = fct_reorder(mod, 
    desc(meanRMSE_test)))

write.csv(mean.xval_axcpt_rt_k_loo.iter, "tables/xval_axcpt_rt_k_loo.mean.csv", row.names = F)

ACC models

all.acc.0 <- glmer(mod.acc.0, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.0.1 <- glmer(mod.acc.0.1, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.1 <- glmer(mod.acc.1, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.1.1 <- glmer(mod.acc.1.1, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.2 <- glmer(mod.acc.2, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.3 <- glmer(mod.acc.3, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

kable(anova(all.acc.0, all.acc.1, all.acc.2, all.acc.3))

Df		AIC	BIC	logLik de	viance	Chisq Ch	i Df Pr	(>Chisq)
all.acc.0	4	13487.86	13519.91	-6739.931	13479.86	NA	NA	NA
all.acc.1	10	13490.65	13570.76	-6735.323	13470.65	9.215386	6	0.1618223
all.acc.2	16	13496.08	13624.26	-6732.038	13464.08	6.570276	6	0.3624199
all.acc.3	25	13495.63	13695.92	-6722.816	13445.63	18.444541	9	0.0303527

kable(anova(all.acc.0, all.acc.1.1, all.acc.2, all.acc.3))

Df		AIC	BIC	logLik de	viance	Chisq Ch	i Df Pr	(>Chisq)
all.acc.0	4	13487.86	13519.91	-6739.931	13479.86	NA	NA	NA
all.acc.1.1	10	13494.19	13574.31	-6737.094	13474.19	5.673929	6	0.4606892
all.acc.2	16	13496.08	13624.26	-6732.038	13464.08	10.111734	6	0.1200248
all.acc.3	25	13495.63	13695.92	-6722.816	13445.63	18.444541	9	0.0303527

kable(anova(all.acc.0, all.acc.0.1))

Df		AIC	BIC	logLik de	viance	Chisq Ch	i Df Pr	(>Chisq)
all.acc.0	4	13487.86	13519.91	-6739.931	13479.86	NA	NA	NA
all.acc.0.1	7	13483.66	13539.74	-6734.830	13469.66	10.20268	3	0.0169196

summ(all.acc.0.1, confint = T, digits = 3)

Observations	22286
Dependent variable	acc
Type	Mixed effects generalized linear model
Family	binomial
Link	logit

AIC	13483.659
BIC	13539.741
Pseudo-R² (fixed effects)	0.299
Pseudo-R² (total)	0.345

Fixed Effects
	Est.	2.5%	97.5%	z val.	p
(Intercept)	1.994	1.911	2.077	46.969	0.000
ConditionAY	-0.332	-0.423	-0.242	-7.162	0.000
ConditionBY	2.430	2.219	2.641	22.597	0.000
L2_exposure	-0.045	-0.109	0.020	-1.353	0.176
aoa	0.015	-0.051	0.081	0.442	0.659
L2_exposure:aoa	0.089	0.027	0.150	2.816	0.005

Random Effects
Group	Parameter	Std. Dev.
participant	(Intercept)	0.484

Grouping Variables
Group	# groups	ICC
participant	455	0.067

summ(all.acc.0, confint = T, digits = 3)

Observations	22286
Dependent variable	acc
Type	Mixed effects generalized linear model
Family	binomial
Link	logit

AIC	13487.862
BIC	13519.909
Pseudo-R² (fixed effects)	0.296
Pseudo-R² (total)	0.345

Fixed Effects
	Est.	2.5%	97.5%	z val.	p
(Intercept)	2.002	1.918	2.085	46.884	0.000
ConditionAY	-0.333	-0.424	-0.242	-7.171	0.000
ConditionBY	2.430	2.219	2.640	22.616	0.000

Random Effects
Group	Parameter	Std. Dev.
participant	(Intercept)	0.495

Grouping Variables
Group	# groups	ICC
participant	455	0.069

anova(all.acc.0)

## Analysis of Variance Table
##           Df Sum Sq Mean Sq F value
## Condition  2 643.32  321.66  321.66

ef <- as.data.frame(Effect(c("L2_exposure", "aoa"), all.acc.0.1, xlevels = list(aoa = c(-2, 
    -1, 0, 1, 2)), confidence.level = 0.686))
ef$aoa <- factor(ef$aoa)
f_acc <- ggplot(ef, aes(x = L2_exposure, y = fit, ymin = lower, ymax = upper, group = aoa)) + 
    geom_line(aes(colour = aoa), size = 1) + geom_ribbon(fill = "grey", alpha = 0.3) + 
    coord_cartesian(ylim = c(0.5, 1)) + ylab("Model-estimated\naccuracy") + xlab("Scaled L2 exposure") + 
    scale_colour_brewer(palette = "Dark2") + theme_bw(base_size = 12) + theme(panel.border = element_blank(), 
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

ggsave("figures/traditional_acc_plot.png", plot = f_acc, width = 12, height = 8, 
    units = "cm")

LME regression analyses (Analysis 1)

Jason W. Gullifer & Debra Titone

2019-07-24

Preliminaries

Loading exec control data

Create folds

Entropy PCA and scaling of variables

Join LHQ and AXCPT

Jobs data

Experiment version

Data plots

Summary stats table

Group level

RT: individual differences

RT: individual differences: raw

ACC: individual differences

Set loo fold

Traditional AXCPT analysis

RT models

RT models cross validated

ACC models